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Abstract. This work is concerned with the derivation of a robust a posteriori error estimator 
for a discontinuous Galerkin method discretisation of linear non-stationary convection-diffusion ini- 
tial/boundary value problems and with the implementation of a corresponding adaptive algorithm. 
More specifically, we derive a posteriori bounds for the error in the L 2 (H 1 )-type norm for an interior 
penalty discontinuous Galerkin (dG) discretisation in space and a backward Euler discretisation in 
time. An important feature of the estimator is robustness with respect to the Peclet number of 
the problem which is verified in practice by a series of numerical experiments. Finally, an adaptive 
algorithm is proposed utilising the error estimator. Optimal rate of convergence of the adaptive 
algorithm is observed in a number of test problems, discontinuous Galerkin methods; a posteriori 
error estimators; time-dependent convection-diffusion equation. 

1. Introduction. The interaction between convection and diffusion modelled by 
initial/boundary value problems involving partial differential equations (PDEs) poses 
a number of challenges in the context of their numerical approximation. Indeed, 
stationary convection-dominated convection-diffusion type problems admit analyti- 
cal solutions of a multiscale nature that can contain steep gradients, usually termed 
boundary or interior layers, depending upon their location in the computational do- 
main. The accurate and efficient numerical resolution of such steep layers is a chal- 
lenge as their exact location cannot be, in general, known a priori. In the special cases 
where the location of boundary or interior layers is known, structured grids have been 
successfully employed ([20 ). Ultimately, in non-stationary convection-diffusion equa- 
tions, the nature of the solution (including layers) may vary throughout the domain as 
time progresses. This renders the use of adaptive algorithms an attractive proposition 
for the accurate and efficient numerical approximation of this class of problems. As 
adaptive algorithms are usually based on suitable a posteriori error estimators, the 
formulation of estimators that can robustly estimate both the temporal and spatial 
nature of the error are of particular interest. 

A posteriori error estimation for stationary linear equations is now relatively well 
understood; for pure diffusion problems there is a huge array of estimators available 
for a wide variety of different types of finite element discretisations ([301 [I]) and for dG 
methods in particular ( |16l l^lfT^]). For stationary convection-diffusion equations, the 
quest for robust a posteriori error estimators (in the sense that they are independent 
of the Peclet number of the problem) has seen recent advancements in various contexts 

(nana in [in [2i OS]). 

A posteriori error estimators for non-stationary linear convection-diffusion equa- 
tions are also available for various discretisations ([HI O EH HOl 12 El [26l H21 HB]). A 
posteriori error estimators for (spatial) dG methods for non-stationary pure diffusion 
problems can be found in [H [TTJ [35] . 

This work is concerned with the derivation and implementation of an a posteriori 
estimator in the L 2 (iJ 1 )-norm for the backward Euler discretisation in time (zero-th 
order dG method) and the interior penalty dG discretisation in space of the non- 
stationary linear convection-diffusion equation with variable coefficients. The deriva- 
tion of the a posteriori bound utilises the elliptic reconstruction technique ( |19l 1 1 3j ) 
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which allows the use of the robust elliptic error estimator from |23j . The temporal 
residual is also treated so as to ensure that robustness with respect to the Peclet 
number of the problem is maintained. The resulting a posteriori error estimator is 
shown numerically to be reliable in both time and space. The a posteriori estimator 
derived below can be viewed as the analogue of the one presented in [28] when the 
spatial discretisation is a dG method. 

The a posteriori estimator is used to drive an adaptive algorithm which is numer- 
ically assessed in two ways. First, the savings with respect to the degrees of freedom 
by using space-time adaptivity are shown in a series of test problems. Second, the 
rates of convergence of the space-time adaptive algorithm are computed. The adap- 
tive algorithm appears to convergence optimally in both space and time when the 
spatial polynomial degree, p, satisfies p > 2. Somewhat surprisingly, for p = 1, the a 
posteriori estimator appears to converge at a suboptimal rate when the mesh changes 
between time steps. The nature of this degeneracy for p = 1 remains unclear to 
the authors: it is in fact possible that the method itself is suboptimal in time when 
mesh change is allowed. It is also unclear if this is a special phenomenon of the dis- 
cretisation presented in this work or a more general phenomenon extending to other 
discretisations of parabolic problems. This phenomenon will be investigated further 
elsewhere. 

The remainder of this work is structured as follows. The function space setting 
and the model problem are given in Section [2] Section [3] describes the discretisation 
of the problem. In Section [4] we state an a posteriori error bound for the stationary 
elliptic problem. Sections [5] and [6] contain the derivation of the a posteriori bounds 
for the non-stationary convection-diffusion model problem in the semi-discrete and 
the fully-discrete settings, respectively. Section [7] describes a space-time adaptive 
algorithm driven by the a posteriori estimators derived in the previous section. A 
series of numerical experiments are presented in Section [H] and the final conclusions 
are drawn in Section [9] 



2. Model problem. Let the computational domain f2 C R 2 be a bounded Lip- 
schitz polygon with boundary dtt. 

We denote the standard L 2 -inner product on Q, by (•, •) and the standard L 2 -norm 
on O by || • ||. For 1 < p < +oo, we define the spaces L p (0,T; X) (where X is a real 
Banach space with norm || • ||x) that consist of all measurable functions v : [0, T] — > X 
for which: 

\\ v \\lp(0.T;X) 
\\ V \\l^(0.T;X) 



( J \\v(t)\\ p x dtj P < +oo, for 1 < p < + 

ssup |jv(i)||x < +oo, for p = +oo. 
:t<T 



(2.1) 



: ess 

0<i<f 



We also define H^-^^X) := {u G L 2 (0,T;X) : u t G L 2 (0,T;X)}. Finally, we 
denote by C(0,T;X) and C o,1 (0, T; X), respectively, the spaces of continuous and 
Lipschitz-continuous functions v : [0,T] — > X such that: 



IMIc(0,T;X) = m t axJKt)|U < °°> 



dv, , P- 2 ) 



MIco.^ctj) = max|||i;||c(o,r;X) 5 II^||c(o,t,x)} < 00 



We consider the model problem of finding u : £1 — > K such that: 



— - eAu + a.-Vu + bu = f infix(0,T], 

u = on aOx(0,T], ( 2 - 3 ^ 
u(-,0) = uq in O. 

We make the following assumptions: u a G L 2 (VL), f G C(0,T; L 2 (&,)), e G (0,1], 
a G C(0,T; W 1,0O (f2)) 2 , and 6 G C(0, T; L°°(0)). 

It is assumed that the length of a = a(x, t) and the area of f2 are of order one, 
possibly up to rescaling, so that e _1 can be taken to be the Peclet number of the 
problem. For simplicity, we assume that there are constants j3 > and c* > such 
that: 

fo-^V-a>/? a.e. in fi x [0,T], || - V • a + 6|| C (o,T,L-(n)) < (2-4) 



The weak form of Q then reads: find u G L 2 (0, T; H^(fi)) n H^O.T; L 2 (Q)) such 
that for each £ G (0, T] we have 

f^vdx + f(sVu ■ Vv + a • Vuv + buv) dx = f fvdx Vv € Hq(Q). (2.5) 

Under the regularity assumptions above, we have that u G C(0, T; L 2 (S1)). 

3. Discontinuous Galerkin method. Let the mesh £ = {K} be a shape- 
regular subdivision of Q, with if denoting a generic element. We assume that the 
subdivision £ is constructed via affine mappings F R '■ K — > K with non-singular 
Jacobian where K is the reference triangle or the reference square. The mesh is 
allowed to contain a uniformly fixed number of regular hanging nodes per edge. We 
define the finite element space 

V h = V h (0 := {v E L 2 (il) : v\ K o F K G V P (K),K G C}, (3.1) 

where V P (K) is the space of polynomials of total degree p if K is the reference triangle, 
or the space of polynomials of degree p in each variable if K is the reference square. 

Denote by £{Q the set of all edges in the triangulation ( and £ mt (Q) ^ £(C) the 
subset of all interior edges. We also denote the diameter of an element K by h R and 
the length of an edge E by He- The outward unit normal on the boundary dK of an 
element K is denoted by n^. Given an edge E G £ mt (Q shared by two elements K 
and K', a vector field v G [H 1 ' 2 ^)] 2 and a scalar field v G H^ 2 (n), we define jumps 
and averages of v and v across E by: 



M = \( v \k + v \k>), M =v|# • n K + v\ R , ■ n K >, 

W = \( v \k +v\r>), M =v\ R n K +v\ R ,n K ,. 



(3.2) 



If 15 C dft, we set {v} = v, [v] = v • n, {v} = v and [v] — vn, with n denoting the 
outward unit normal to the boundary dfl. 

We define the inflow and outflow parts of the boundary dft at the time t, respec- 
tively, by: 

dn\ n = {x G dn | a(x,t)-n(x) < 0}, drt out = {x e dn \ &(x,t)-n(x) > 0}. (3.3) 



The inflow and outflow parts of an element K at time t are similarly denned as: 



dK\ n = {x £ dK | a(ar, t) ■ n K (x) < 0}, dK l out = {x£dK\ a(x, t) ■ n K (x) > 0}. 

(3.4) 



The semi-discrete discontinuous Galerkin approximation to (2.5) then reads as 
follows. For t = 0, set Uh(0) £ Vh to be a projection of uo onto Vh- For t £ (0, T], 
find u h £ C°^(0,T; V h ) such that 

du 

(-gj-,v h ) + B(t;u h ,v h ) + K h (u h ,Vh) = (f,Vh) Vvh£Vh, (3.5) 

with 

B(t;w,v) =B(£;t;w,v) 



:= / (eVu> — aw) ■ Vv + (b — V • a)wv dx + ~ \w]-\v]ds 

Ik h-E lw 

KeC K Ee£(C) 

+ /J / a ■ iiR-tin) ds + / a ■ n K w(v\^ — v\^,) ds \ , 

^(lo, «) =i^(C; to, «) := - 51 / i £ ^ w } ■ M + { eVu i ■ N ds > 

Ee£(C) E 

(3.6) 

where K OK' — E C dKl ut \dfl. In standard fashion, the penalty parameter 7 > is 
chosen large enough so that the operator B+Kh is coercive. Moreover, for simplicity of 
the presentation, we assume 7 > 1 so that the constants in the subsequent discussion 
are independent of it. 

We note that the bilinear form Kh is not well-defined for arguments in H^(fl), 
but the bilinear form B is and for u, v £ H^fl) and t £ (0, T], we have 

B(t;u,v)= / (eVu • Vw + a • Vuv + buv) dx. (3.7) 
Jn 



In light of this, the weak problem (2.5 1 can be rewritten for each t £ (0, T] as 



(ill 

(—,v)+B(t;u,v) = (f,v) Vv£Hi({l). (3.8) 



We shall also consider a full discretisation of problem (2.5) by using a backward 
Euler method to approximate the time derivative. 

To this end, consider a subdivision of [0, T] into time intervals of lengths t\ , T2, t„ 
such that Y^=i T j = T for some n > 1 and set f° = and t k := Ylj=i T k- Denote 
an initial triangulation by (°. We associate to each time step k > a triangulation 
£ fe which is assumed to have been obtained from C by locally refining and coars- 
ening £ . This restriction upon mesh change is made to avoid degradation of the 
finite element solution, cf. [9]. To each mesh ( k , we assign the finite element space 



V£ = V h (C k ) given by pT) 



We set B k (t; •, •) := B(( k ; t; •, •) and K k (-, •) := K h (( k ; ; ■). Moreover, when the 
last two arguments of B are in Hq(Q,), we simply denote the bilinear form by B. We 
also denote f(-,t k ) — f k , a(.,t k ) — a k , and b(.,t k ) — b k for brevity. 



The fully-discrete dG method then reads as follows. Set it° to be a projection of 
u onto V,°. For k = 0,...,n - 1, find u k h +l G V h fc+1 such that 



Tk+l 



V^ +1 G V h k+1 . 

(3.9) 

We shall take u h to be the orthogonal L 2 -projection of uq onto V°, although other 
projections onto V® can also be used. 

4. Error bounds for the dG method for the stationary problem. To 

analyse the error, we introduce the following quantities: 

MllcH E( £ ll v *w+«W)) + E SlMlliW'" 



i«kc:=ff Jn ":;^ )'+ e K+-)iiMiii^) 

VVt-eirjCnjuo} IIMII / Ee£(Q v e 7 



1/2 



(4-1) 

We note that ||| • |||^ and | • \ a,( define norms on Hq(£1) + Vh- When no confusion is 
likely to occur we shall suppress the dependence on the mesh £ and write 1 1 1 • 1 1 1 and 
Ha- 

Throughout this work, the symbols < and > are used to denote inequalities true 
up to a positive constant independent of e, f3, u, and Uh- 

Let t e (0, T] be fixed. It is easy to see that the bilinear form B(t; •, •) is coercive 
on Hq(CI), viz., 

B(t;v,v)>\\\v\\\ 2 , (4.2) 

for all v £ Hq(£1), and is continuous in the following sense: 

B(t;w,v)<(\\\w\\\ + \w\ A )\\\v\\\, (4.3) 

for all w € Hqity+Vh, and v E i?Q (f2). Moreover, the discrete bilinear form is coercive 
for Vh in Vh with respect to the ||| • ||| norm, viz., 

B{t;v h ,v h ) + K h (v h ,v h ) > IIKIH 2 (4.4) 

Next, we introduce the following notation which will be used to define the a posteriori 
estimators: 

olk ■= mm(h K e~^ , /3 - ^), ole '■= min(/i £ e~5 ; /3~3), ar = min(e~5 j p-k ). (4.5) 

An a posteriori estimator for the stationary problem inspired by 23J will be utilised in 
our analysis. More specifically, we have the following result whose proof is completely 
analogous to the one of Theorem 3.2 in [53] and is therefore omitted for brevity. 
Theorem 4.1. For a given t € (0,T], let u s be such that 

B(t;u s ,v) = (f,v) Vw 

and consider u s h € Vh such that 

B(t;u s h ,v h ) + K h (u s h ,v h ) = (f,v h ) \/v h e V h . 



Then the following a posteriori bound holds: 



<\\\ + K - <\a) 2 < E a KWf + eA < - a • V< - bu s h \\ 2 L2{K) 

+ e ^^m<}\\U E) + E {ir + P hE + h T 



l h\\\L 2 (E)- 



5. An a posteriori bound for the semi-discrete method. To highlight the 
main ideas, we begin by deriving an a posteriori bound for the semi-discrete problem. 

Definition 5.1. For each t e [0,T], we define the elliptic reconstruction w £ 
Hq(Q) to be the (unique) solution of the problem 

B(t;w,v) = (f-^,v) VveHfcSl). 



Remark The dG discretisation of Definition 5.1 is to find a function Wh G C ' 1 (0, T, Vh) 
such that for each t £ (0,T] we have 



B(t; w h ,v h ) + K h (w h ,v h ) = (/ - 



duh 
dt 



,vh) Vv h eVh 



This, in conjunction with (3.5) and (4.4), implies that Wh = u^. Thus, |||u> — Uh 



\w — uh\a can be estimated using Theorem 4.1 
We decompose the error as follows: 



e = u — Uh = p + 9, with p := u — w and 9 := w — Uh- 



(5.1) 



The dG solution, Uh, admits a decomposition into a conforming part Uh. c G i^Q (r2)n V/,. 
and a non-conforming part u^^d € Vh with Uh = Uh. c + Uh.d, such that: 



\ u h,d\ 



+ \u h4 ?A< E ff I + + — )IIK]II 



\h E 



h E 

e 



1-^11 Z 1^ h E \\\—\\\ L , {Ey 



(5.2) 



We refer to [23] for proof of these estimates which are based on respective constructions 
by [TB] and [25]. We further define e c :— u — Uh, c and 9 C := w — Uh. c - 
Lemma 5.2. For each t € (0, T] and for all v e Hq(CI) we have 

de 

( — ,v) + B(t;p,v) = Q. 



Proof. This follows directly from Definition 5.1 and (3.8). □ 



We define the initial condition estimator, fji, by 
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\u - u h 



E MIK 

Ee£(C) 
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i 2 (B)' 



(5.3) 



while the spatial estimator, fjs, is given by 



V 2 s 



J Vsi dt + min { ( J Vs 2 dt ) >4 J Vs 2 dt }> 



(5.4) 



where 



VSt = E a K\\f- -kt +£&Uh -a- Vu h ~ bu h \\ 2 L2{K) + E £ia E\\[Vu h }\\ 2 L2{E) 



Kec 



E { E f L +PhE+ h f-)\\[u h ]\\h (EY 
Ee£(C) E 



is the residual term and 



= E * 



du h 2 



(5.5) 
(5.6) 



is the rate of change of non-conformity. 

Theorem 5.3. The error e of the semi- discrete dG method (3.51 satisfies the 
bound 



[ l||e||| 2 <ft<^+^I. 
Jo 



Proof. Choosing v — e c in Lemma [5 . 2 1 gives 

e c ) + B(t; e c , e c ) = (—^' e c) + B (t\ &c, e c ). 



Using (4.2 1, (4.3), and Young's inequality we arrive to 



dt 



|e c || 2 ) + |||e c 



< 



+ |0cU) 2 + ||^| 



(5.7) 



(5.8) 



Let E c := ||e c (T )|| with T E [0,T] such that ||e c (T )|| = max < t < T ||e c (i)|| then 
integrating on [0, T] gives 



2 dt< ||e c (o)|| 2 + / (\\\e c \\\ + \e c \ A ) 2 dt + E c 



at 



I dt. (5.9) 



Integrating (5.8 1 on [0, To] and using Young's inequality once again gives 

(£ c ) 2 <IM0)|| 2 + / T (|||0 c ||| + |0 cU )2 d t + (^ T ||^|| dt ) 2 . 







(5.10) 



Using (5.10) to bound E c on the right-hand side of (5.9) yields 

pT 

\ 2 dt<. 



+ I' (lll^lll + \e c \ A ) 2 dt+ ( fw^Wdtf. (5.ii) 



Going back to (5.8), using the Poincare-Friedrichs inequality on ||e c || and working as 
above gives 



I e r 1 1 1 2 dt < lie, 



+ I (\\\9 c \\\ + \9 c \ A ydt + a 2 T 
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I du h ,d I 
1 0i 1 



(5.12) 



Using the triangle inequality, we have 
|||e||| < |||e c ||| + Wu/j.dlll and |||0 £ 



+ \O c \a < \\\d\\\ + \0\A + \\\u h , d \\\ + \u hid \ A . 

(5.13) 



Combining ( |5.11[ ), ( |5.12| ), and ( |5.13[ ) yields 



|e||| 2 di < He, 



o 



{( 





T n du h} , 



|0|a) dt+ I (\\\u h ,d\ 
>o 



\uh,d\ 2 A) dt 



\'H) -<-c, / fT ||^|| 2 ^}- 



(5.14) 



The proof then follows directly from Theorem 4.1 and (5.2). 
□ 

6. An a posteriori bound for the fully-discrete method. We now con- 
tinue by applying to the fully-discrete setting the general framework presented in the 
previous section. 

Definition 6.1. We define the elliptic reconstruction w k+1 £ Hq(Q) to be the 
unique solution of the elliptic problem 



B{t 



V) = (f 



Tk+1 



At each time step k, we decompose the dG solution u k into a conforming part 
u h,c e #o(^) n Vh an d a non-conforming part u\ d £ V k such that u 1 ^ = u k hc + u\ d . 
Given t £ (t k , t we (re) define Uh(t) to be the linear interpolant with respect to t 
of the values and u k+1 , viz., 



u h {t) := l k {t)u k h + l k+1 (t)u k+1 , 



(6.1) 



where {lk,lk+i} denotes the standard linear Lagrange interpolation basis defined on 
the interval [t k ,t k+1 \. We define Uh. :C (t) and Uh,d{t) analogously. We can then decom- 
pose the error e = u — = e c — Uh,d where e c = u — it^.c- It will also be useful to 



define 9 k 



Lemma 6.2. Given t £ (t k ,t k 



k +fc+ll 



have 



,de 
di' 



(£» + B(t; e c , v) = B(t k+1 ;w k+1 ,v) - B(t, u h , c , «) + (/- f k+ \ v) Vt> £ H%(Sl). 



Proof. This follows from Definition 6.1 and (3.8 1 and by rearranging the resulting 
equation. □ 

Before proving the a posteriori bounds for the fully-discrete method, we introduce the 
error estimators. We begin by defining the initial condition estimator, rjj, by 



Vl = IK-u"|| 2 + h E\\[Uh\\\h(E), 
EG£(C°) 



(6.2) 



while the spatial estimator, r/s, is given by 



3=0 



3=0 



n-1 r t i+1 



n-1 ~ti+ 1 ^2 -t-.- 



where 



i+i j 
u h — u 



L 2 (K) 



Be£ i " t (C 3 uC 3+1 ) Ee£(Cu£'+ 1 ) 
£e£(C j uc j '+ 1 ) £ 

and 

Be£(C J 'uC J+1 ) 3+ 
+ J] ^ X |lfeWaK +1 - <) + (a J+1 - a)< +1 ]||i 2(E) . 

Ee£i«t(£JUC' +1 ) 

The final term is the time (or temporal) estimator, tjt, given by 

n— 1 /-n— 1 / .t J + 1 \2 n-l / ,t J + 1 



(6.4) 



(6.5) 



^n—x x ft \* ft 

2 E t j+i^ti j+i + min { ( E / 'few+i dt ) '^E / . ^ • 

(6.6) 



where 



^■+1= E e||V(< +1 -<)||| 2(/c) , (6.7) 



and 



4 2J -+i = E iM')( a • - <) + & K +1 - O) + f f j+1 

A'e^uC^ 1 (6.8) 
+ (a^ 1 - a) • V< +1 + b)u{ +1 \\l 2{K) . 

Theorem 6.3. The error e of the fully-discrete method satisfies the bound 

n-l „ t i+i 

E / IINII 3 *^^ +Vs + Vt- (6-9) 

3=0 ^ 

Proof From Lemma |6.2| we have 



(6.10) 

which upon straightforward manipulation and setting v = e c gives 

~(| |e c || 2 ) + e c , e c ) =(^, e c ) + B(t k+1 ; 9 k c + \ e c ) + l k (t)B(t; - < c , e c ) 

+ B^ufc 1 , e c ) - B(t; ^+ c \ e c ) + (/ " /** S e c)- 

(6.11) 
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Denoting by l\e c the L 2 -projection of e c onto H^(Q) H V h {( k U C k+1 ), using ( |42| ), 
(4.3), the Cauchy-Schwarz inequality, Young's inequality and the triangle inequality 



yields 

d , 



•(|||^ +1 ||| + |^ +1 U) 2 + IIK rf ||| 2 + K rf |^ 
t 1 + t 2 + t 3 + t a + t 5 , 



where 



Ti := ^ / (/ - / fe+1 + (a fe+1 - a) • V< +1 + (fe fe+1 - b)u k h +l 



+ Z fc (t)(a • V(< +1 - «*) + 6« +i - uf)))e c dx, 



k+1 



(6.12) 



(6.13) 



T 2 := hit) J2 I eV K +1 - u t) ■ Ve c dx, 



(6.14) 



To 



A'fc 



[(a fc+1 -a)< +i ]( ec -4 fc e c )d S , 



,fe+i 



(6.15) 



(is, 



and 



T B := 



5 E 



A'ec fc 



fe (t)a« +1 - 4) + (a fc+1 - a)< +1 ]4 te e C d S . (6.17) 



(6.16) 



Using the discrete and integral versions of the Cauchy-Schwarz inequality, we imme- 
diately have 



T\ < ?7T 2 ,fc+i||e c ||, T 2 < r?Ti,fe+i||l e c 
Further, using the approximation properties of it we have 

T 3 + Ti < r?Si,fc+i|||e c |||- 



(6.18) 



(6.19) 



Next, using the standard inverse estimate H-^cIIi^e) < ^^ll-^clll^A')' f° r E C 
dK, and the stability of the L 2 -projection ||/*e c || < ||e c ||, we deduce 



T 5 < VS 2 ,k+i\\e c 



(6.20) 



The proof then follows from Theorem 4.1 and (5.2) and by employing a bounding 



strategy identical to that used in Theorem |5.3| □ 



RemarkAn alternative a posteriori bound can be show, using directly (4.3) on 



(6.11). The resulting bound will contain computable terms in the | • l^-norm, which 
can, in turn, be estimated by a more easily localisable term of the form e _1 / 2 || • ||i,»(nv 
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Numerical experiments (omitted here for brevity) show that this alternative estimator 
produces very similar results to the ones based on Theorem |6.3| 

RemarkAlthough efficiency bounds are not presented in this work, we refer to [2 8) 
for the proof of efficiency bounds for the a posteriori estimator presented therein for 
the SUPG spatial setting which should be extendable to the estimator presented in 
Theorem |6.3| given the similarity of some of the terms presented. 

RcmarkThe use of elliptic reconstruction is not essential for the proof of Theorem 
|5.3| and Theorem |6.3| It is possible to derive the residual based a posteriori bounds 
directly albeit at the cost of a lengthier calculation. The advantage of using the 
elliptic reconstruction in the proof lies in the fact that it can be easily modified to 
include non-residual based spatial a posteriori estimators. This, in turn, may offer 
improvements in robustness with respect to the Peclet number (cf. [H]). 

7. An adaptive algorithm. The a posteriori bounds presented above will be 
used to drive a space-time adaptive algorithm. A number of adaptive algorithms 
for parabolic problems have been proposed in the literature; see e.g. HI) and the 
references therein. 

Here, we propose a variant of the adaptive algorithms from [51 [22] which appears 
to perform well for our discretisation. The pseudocode is given in Algorithm 1. The 
algorithm is based on using different parts of the a posteriori estimator from Theorem 



(6.3) to drive space-time adaptivity. 

Both mesh refinement and coarsening are driven by the term r]s 1 j+i- The size of 
the elemental contributions to rjs ly j+\ determines the percentages of elements to be 
refined and/or coarsened depending on the two spatial tolerances stola and stolb. 
The nature of the time estimator, r]T, makes it difficult to use as a temporal refinement 
indicator so, to this end we define f]T,j+i given by 

1 f tJ+1 

Vt,j+1 = l r j+lV Tl ,j+l + min {"T,r} J ^ T)T 2 ,j+l dt - C 7 - 1 ) 

The sum of all terms ?7tj+i bound r\T and hence can be used to drive temporal 
refinement subject to a time tolerance ttol on each time interval. 

RemarkThe term rjs 2 j+i is smaller than rjs 1 j+i both from a theoretical and a 
practical standpoint. It is, therefore, not used to drive the spatial mesh modification. 

A difference of Algorithm 1 to (standard) adaptivity algorithms for conforming 
finite elements for parabolic problems is that here a limit is given on how much 
the mesh can change. This is introduced in order to to ensure termination of the 
algorithm. The input value m is a bound on how many times the mesh can change 
per unit time. The value of m can be very large in practice. 

8. Numerical experiments. We shall investigate numerically the presented a 
posteriori bounds and the performance of the adaptive algorithm through an imple- 
mentation based on the deal . II finite element library ( |4J). All the numerical exper- 
iments have been performed using the high performance computing facility ALICE at 
the University of Leicester. 

We denote by the total number of degrees of freedom on the union mesh 
C k U C fe+1 - Hence, the weighted total number of degrees of freedom of the problem is 
given by 

n-i 

Total DoFs := ^ r j+1 Xj. (8.1) 
o=o 
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Algorithm 1 Space-time adaptivity 



1: Input: e, a, 6, /, uq, T, Vl, m, n, C°, 7, initol, ttol, stola, stolb, ref %, coar%. 

2: Calculate v? h . 

3: Set j — 0, Ti, r„ = T/n, time = n, threshold — T/m, counter = 0. 

4: while 77/ > initol do 

5: Refine C° by 10% and coarsen C° by 5%. 
6: Calculate 

7: end while 

8: while izme < T do 
9: Calculate from u{. 
10: while f/T.j+i > ttol do 

11: n <r- n + 1 

12: time ■<— time — r J+ i 

13: T J+3 = T J+2j •••! T n = T n-1 

14: T i+ 2 = Tj+i/2 

15: T i + i «- T, + i/2 

16: time + Tj + i 

17: Calculate from ui. 

18: end while 

19: counter <— counter + r^+i 

20: if r)s 1 ,j+i > stola and counter > threshold then 
21: Create from & by refining by ref %. 

22: end if 

23: if stolb < r)Si,j+i < stola and counter > threshold then 

24: Create from by refining by ref % and coarsening by coar%. 

25: end if 

26: if ?7si,j+i < stolb and counter > threshold then 
27: Create from by coarsening by coar%. 

28: end if 

29: Calculate from u 3 h . 

30: j ir- j + 1 

31: time <— time + Tj +1 

32: if counter < threshold then 

33: counter = 

34: end if 

35: end while 



In all examples presented below, unless otherwise stated, 7 = 10, stolb = stola/5, 
ref% = 6.25% and (° is an 8 x 8 uniform quadrilateral mesh. We further take 
coar% = 10% for Example 1, Example 2 and Example 4 while we set coar% = 30% 
for Example 3. 

8.1. Example 1. Let Q = (0, l) 2 , a = (1, 1) T , b = 0, u Q = 0, T = 10 and select 
the function / so that the exact solution to problem (2.5) is given by 

u(x, y, t) = (1 - e"*) { ' e _ 1/e - 1 { l + y - lj . (8.2) 

The solution exhibits boundary layers at the outflow boundary of the domain of width 
0(e) as well as a temporal boundary layer. 
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Time steps 


Total DoFs 


Est. 


Est. Ratio 


Error 


Err. Ratio 


10 


160 


3.83e-l 




6.45e-2 




20 


640 


1.94e-l 


0.507 


3.07e-2 


0.475 


40 


2560 


9.79e-2 


0.502 


1.47e-2 


0.479 


80 


10240 


4.90e-2 


0.501 


7.176-3 


0.487 


160 


40960 


2.45e-2 


0.500 


3.53e-3 


0.492 


320 


163840 


1.22e-2 


0.500 


1.756-3 


0.496 


640 


655360 


6.14e-3 


0.500 


8.72e-4 


0.497 


1280 


2621440 


3.07e-3 


0.500 


4.35e-4 


0.498 



Table 8.1 



Example 1. Decay of the a posteriori error estimator for p = X, e = 1. 



Time steps 


Total DoFs 


Est. 


Est. Ratio 


Error 


Err. Ratio 


10 


160 


1.66c+l 




8.56e-l 




20 


640 


1.29c+l 


0.773 


1.06c+0 


1.240 


40 


2560 


1.08e+l 


0.838 


1.18c+0 


1.114 


80 


10240 


8.68e+0 


0.802 


1.27e+0 


1.078 


160 


40960 


5.91e+0 


0.680 


1.14e+0 


0.899 


320 


163840 


3.47e+0 


0.587 


7.54e-l 


0.657 


640 


655360 


1.88e+0 


0.542 


4.06e-l 


0.538 


1280 


2621440 


0.98e+0 


0.520 


2.01e-l 


0.496 



Table 8.2 



Example 1. Decay of the a posteriori error estimator for p = 1, e = 10 



We begin by assessing the decay of the estimators on uniform spatial and temporal 
meshes. To this end, the number of time-steps is doubled while the number of spatial 
degrees of freedom is quadrupled by respective uniform refinements. The results for 
the spatial polynomial degree p = 1 are given in Tables |8.1| and |8.2| As expected, 
the error halves (Err. Ratio) for each space-time uniform refinement and the same 
is observed for the a posteriori estimator (Est. Ratio), implying optimality of the a 
posteriori error estimator for this example. This is achieved immediately for e = 1 
and upon sufficient resolution of the boundary layer for e = 10 -2 . The a posteriori 
estimator appears to be of optimal rate under uniform refinement in all the numerical 
examples presented in this work; these results are omitted for brevity. 

We now assess the performance of the space-time adaptive algorithm presented 
in Section [7j we do this for p > 1 as the case p = 1 needs special care and hence will 
be discussed separately. 

We begin by fixing a temporal tolerance that produces enough time steps so that 
the temporal contribution to the error is very small in comparison to the spatial 
contribution. The spatial tolerance is then gradually reduced to compare the rates 
of convergence of the spatial estimator with that of the actual total error. Results 
for e = 1 and e = 10~ 2 and for polynomial degrees p — 2 and p = 3 are shown in 



Figure 8.1 here as in all the following convergence plots an unmarked line indicates 



the optimal order of convergence. Optimal rates of convergence are observed with 
respect to the total degrees of freedom for both the estimator and the error. Further, 
as shown in the bottom plots of Figure |8.1| the effectivity indices are bounded and 
remain between 6 and 11 for the two values of e; these are directly comparable to 
those observed in [23J for the stationary problem. 
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Sparjal EstmaloriErrar for n=10 and p=2 



I... 



10 




10 



10 



i 10 



! 10 



10 



Spatial EstimalorfErrot (or i:=10 and p-l 




10 10' 
Total DoFs 



10 10 
Total DoFs 



Spatial Estimator/Error for ;:=1 and p=2 




Spatial Estamatar/Error for :=1 and p=5 




Spatial Effeavily Indices for si= 1 



10* 10* 
Total DoFS 




Fig. 8.1. Example 1: Spatial rates of convergence and effectivity indices for e = l,e = 10 2 
and p = 2, 3 



In order to study the temporal rates of convergence of the adaptive algorithm 
any boundary layers must be fully resolved so that the spatial error is dominated by 
the temporal one. To this end, we set 7 = 100 and use p = 10 on a sufficiently fine 
fixed spatial mesh. We monitor the algorithm's behaviour upon gradually reducing 



the temporal tolerance. The results are given in Figure 8.2 Optimal order is observed 



with respect to both the time estimator and the error. The effectivity indices appear 
to remain bounded and converge to a smaller value for e = 10~ 2 than for e = 1 which 
is evidence of robustness with respect to the Peclet number. 

The presence of a temporal boundary layer in the solution motivates a comparison 
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Temporal Estimator/Error for i:=1 




Temporal Estimator/Error for b"1 



10 10 

Time Steps 




10" 10 10 

Time Steps 



Temporal Effecivitv Indices for i:=10"^ 



Temporal Effecivity Indices for e"1 




10 10 

Time Steps 




Fig. 8.2. Example 1: Temporal rates of convergence and effectivity indices for e = 1, 10 



Time Estimator for ;-.= 1 0" Adaptive vs Uniform 



Time Estimator for G _ 1 Adaptive vs Uniform 




& Adaptive 
V Uniform 




10" 10 
Time Steps 



Fig. 8.3. Example 1: Adaptive vs uniform comparison for e = 1, 10 



between adaptive and uniform time-stepping. We do this by using p = 3 and a suffi- 
ciently resolved spatial mesh and by reducing the temporal tolerance in an adaptive 



fashion and uniformly, respectively. The results are given in the plots in Figure 8.3 
the savings achieved by using adaptive time-stepping are evident. 

Let us now discuss the behaviour of the space-time a posteriori estimator in the 
case p = 1. Recall that the estimator has already been proven to be of optimal rate in 
both space and time under uniform refinement, cf. Tables [8~T| and [872] In the adaptive 
setting, however, where the mesh can change between time steps, we observe that the 
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Fig. 8.4. Example 1: Spatial and temporal rates of convergence for e = 10 2 and p = 1. 



rate of convergence of the time estimator for p = 1 may be slower than the expected 
optimal rate. This is indeed the case for this test problem, as shown in the right plot 
of Figure |&4| by comparison with the results obtained with p = 2 and p = 3. On the 
other hand, the spatial estimator retains optimal rate, as shown in the left plot of 
Figure |8.4| 

Notice that the method itself could be suboptimal under mesh change for p = 1; 
this is something we have not been able to test and indeed the nature of this degen- 
eracy remains unclear to the authors. It is also unclear if this is, or not, a special 
feature of the backward-Euler-dG method considered in this work. In the literature, 
numerical experiments including rates of convergence of space-time adaptive algo- 
rithms for parabolic problems are very rare, especially for non-conforming methods. 
Therefore, it is not possible at this point to assess the generality of this phenomenon. 
We mention in passing the recent work [5] where degeneracy due to mesh modification 
(including mesh refinement) has been observed for Crank-Nicolson FEM discretisa- 
tions. As shown above, the degeneracy of the temporal rate of convergence when 
spatial mesh-change is present disappears for p > 1. Hence, for the remainder of this 
work, all the numerical experiments with space-time adaptivity are performed with 
p > 1. 

8.2. Example 2. We set = (-1, l) 2 , a = (1, 1) T , 6=1,/ = sin(5t)xy, u Q = 
and T — 2tt. The solution exhibits layers of width 0(e) in the proximity of the outflow 
boundary and is oscillatory in time. The sharpness of the boundary layers depend on 
time, thus making this a good test of the ability of the algorithm to add and remove 
degrees of freedom. 

As in Example 1, we begin by fixing a temporal tolerance while decreasing the 
spatial tolerance to observe the rates of convergence for the space estimator. We 
then set p = 3 with a spatial tolerance small enough to resolve any boundary layers, 
while reducing the temporal tolerance to observe the rates of the time estimator. The 
results are displayed in Figure |8.5| Optimal rates of convergence are observed for 
both the space and the time estimators. 

To assess the mesh change driven by the adaptive algorithm we also plot the 
individual degrees of freedom on each mesh against time for a given spatial tolerance 
and temporal tolerance. The results are given in Figure |8.6| We observe that the 
adaptive algorithm is adding and removing degrees of freedom at a rate that is in 
accordance with the oscillating nature of the solution driven by the sinusoidal forcing 

16 



Spatial Estimator for i:=10 




10 
10 

I 10 " 1 

E 

uj 10 

n 

UJ , 

.=■10° 
10 
10 



Spatial Estimator for e»1 



■P=2 
P=3 




Time Estimator for c-1 0" 



10 

Total DoFs 



Time Estimator for c=1 




Fig. 8.5. Example 2: Spatial and temporal rates of convergence for s = 1, 10 2 . 

DoFs vs Time for l-1 



DoFs vs Time for .=10* 




1 B000 





Fig. 8.6. Example 2: Degrees of freedom on each time step plotted against time for e = 1, 10 



function /. 

8.3. Example 3. Let £1 = (-1, l) 2 , T = 100, a = (y, —x) T , b = 0, / = 0, and 
u = e^ 64 ^- - 5 ) e -64(y-0.5) ppjE convects the initial two dimensional Gaussian 

profile along the circular wind while diffusing it at a rate depending upon e. 

Setting p = 3 and fixing a spatial tolerance and an initial tolerance, we reduce 
the temporal tolerance to observe the temporal rates of the problem. Additionally, 
for e = 1, we compare adaptive time-stepping and uniform time-stepping. The results 



are given in Figure 8.7 Some meshes at various time steps produced by the algorithm 
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Time Estimator t: Comparisons 




Time Estimator for i:=1 Adaptive vs Uniform 




10 10 
Time Steps 

Fig. 8.7. Example 3: Convergence of the time estimator for e = 1,10 _2 ,10 — 3 and comparison 
of adaptive vs. uniform time steps. 



Grid at t=0 



Grid at t=2 



Grid at t=4 



Grid al t=6 



Fig. 8.8. Example 3: grid snapshots. 



for e = 10 3 are shown in Figure 

The smaller the value of e, the longer it takes the time estimator to reach optimal 
order which is to be expected theoretically. Furthermore, large savings are observed 
using adaptive time-stepping over uniform time-stepping when the diffusion is large; 
time steps would be otherwise wasted when the solution is effectively zero. 

8.4. Example 4. Let Q, = (0,1) 2 , a = (sin(t) , cos(t)) T , b = 0, / = 1, u = 
and T = 2n. The nature of the solution is rather uniform in time but has a boundary 
layer of width 0{e) whose location does depend on time. Therefore, this example is 
well suited to test the ability of the algorithm to adapt the grid to the feature of the 
solution as time evolves. 
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Spatial Estimator for c-1 




Spatial Estimator for 



10 10 
Total DoFs 




10 10 
Total DoFs 




Time Estimator for c=1 




Fig. 8.9. Example 4- Spatial and temporal rates of convergence for e = 1, 10 . 

As in previous examples, we fix a temporal tolerance and then reduce the spatial 
tolerance to observe the rates of the space estimator. Again, we also fix p — 3 and 
a spatial tolerance small enough to ensure that all boundary layers are sufficiently 
resolved and then observe the rates of the time estimator. These results are given 
in Figure 8.9 Finally, grids at various times for e = 10~ 2 and p 



Finally, grids at various times for e = 10 -2 and p = 2 are given in 
Figure |8. 10 Optimal spatial and temporal rates of convergence are observed. The 
grids produced for e = 10~ 2 clearly show that the adaptive algorithm is picking up 
the boundary layers as they move around the domain and that unneeded degrees of 
freedom are not retained. 

9. Conclusions. An a posteriori error estimator for the discontinuous Galerkin 
spatial discretisation of a non-stationary linear convection-diffusion equation is pre- 
sented. The derivation of the estimator is based on reconstruction techniques to 
make use of robust a posteriori estimators for elliptic problems already developed in 
the literature. Our numerical examples clearly indicate that the error estimator is 
robust and the respective space-time adaptive algorithm works well for the studied 
problems. The spatial effecitivity indices are in an identical range to those observed 
in (|23j) and the temporal effectivity indices are approximately of the same order as 
those currently seen in the literature for similar problems (|13j) and appear indepen- 
dent of e. Further study to understand the temporal rate of convergence for p = 1 
under mesh-modification as well as extensions to high order time-stepping methods 
will be considered in the future. 
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